%------------------------------------------------------------------------%
%              Comunicaciones digitales y analógicas                     %
%                      Trabajo Práctico N 1                              %
%                                                                        %
%                       Eric Lionel Koplin                               %
%------------------------------------------------------------------------%
close all;clear all;clc;

%% 2) Probabilidad de error
%---------------------------------------------------------------Sim
cant_modulaciones = 4;
cant_constelaciones=4;
for modindex=1:cant_modulaciones    % distintos moduladores
    for i=1:cant_constelaciones     % distintas cantidades de símbolos en el modulador

        M = 2^i;        % cantidad de simbolos
        k = log2(M);    % bits por simbolo
        SNRbdB = 0:10;  % SNR por bit en dB
        nsimbs = 2e4 ;  % cantidad de simbolos que para cada valor de SNRbdB en el Montercarlo

        %----------------------------------- Conversiones e inicializaciones varias
        %
        SNRbveces = 10.^(SNRbdB/10); % SNR por bit en veces
        SNRdB = SNRbdB + 10*log10(k);% SNR en dB
        SNR = 10.^(SNRdB/10);% SNR en veces

        Eavg = 1 ;
        Pes = zeros(1,length(SNRbdB));% probabilidad de error por simbolo simulada
        Peb = zeros(1,length(SNRbdB));% probabilidad de error de bit simulada

        %-----------------------------------------------------Simulacion Montecarlo
        
        switch modindex
            case 1
            %--------------------PAM
            % modulador
            modulator = comm.PAMModulator(M,'SymbolMapping','Gray',...
                                      'NormalizationMethod', 'Average Power',...
                                      'AveragePower',Eavg);
            %demodulador
            demodulator = comm.PAMDemodulator(M,'SymbolMapping','Gray',...
                                          'NormalizationMethod' , 'Average Power' , ...
                                          'AveragePower',Eavg);
            %P(error) teorica :
            % Probabilidad de error de bit teorica de una M−PAM
            % Uso un puntero a una funcion cuyo parametro es SNR (en veces)
            Pet = @(SNR) 2*(1-1/M) * qfunc(sqrt(3*SNR/(M^2-1))) ;                         
            %
            modulator_name=sprintf('%i-PAM',M);
            case 2
            %--------------------PSK
            %Recordar la energia está normalizada, Es=Eavg=1
            % modulador
            modulator = comm.PSKModulator(M,'SymbolMapping','Gray');% Eavg=1
            %demodulador
            demodulator = comm.PSKDemodulator(M,'SymbolMapping','Gray');
            %P(error) teorica :http://www.federica.unina.it/ingegneria/digital-communication/error-probability-for-m-psk/
            %Pet = @(SNR) 2* qfunc(sqrt(2*SNR)*sin(pi/M));                         
            Pet = @(SNR) 1-quadv(@(theta) cos(theta)*sqrt(SNR/pi).*exp(-SNR.*sin(theta)^2),-pi/M,pi/M);                         
            %
            modulator_name=sprintf('%i-PSK',M);
            case 3
            %--------------------QAM
            % modulador
            modulator = comm.RectangularQAMModulator(M,'SymbolMapping','Gray',...
                                      'NormalizationMethod', 'Average Power',...
                                      'AveragePower',Eavg);
            %demodulador
            demodulator = comm.RectangularQAMDemodulator(M,'SymbolMapping','Gray',...
                                          'NormalizationMethod' , 'Average Power' , ...
                                          'AveragePower',Eavg); 
            %P(error) teorica :
            Pet = @(SNR) 1-(1-2*(1-1/sqrt(M)) * qfunc(sqrt(3*SNR/(M-1)))).^2;                         
            %
            modulator_name=sprintf('%i-QAM',M);
            case 4
            %--------------------FSK
            %P(error) teorica :
            Pet = @(SNR) 1/sqrt(2*pi)*quadv(@(x)(1-(1-qfunc(x)).^(M-1))*exp(-(x-sqrt(2*SNR)).^2/2),-1000,1000) ;                         
            %
            modulator_name=sprintf('%i-FSK',M);
        end
        %--------------------AWGN Channel
        %hAWGN = comm.AWGNChannel('NoiseMethod','Signal to noise ratio (Eb/No)',...
        %    'EbNo',SNRbdB(1),'SignalPower',Eavg,'BitsPerSymbol',k);
        %hAWGN = comm.AWGNChannel('NoiseMethod', ...
        %                         'Signal to noise ratio (SNR)','SNR',0);
        %--------------------Create an error rate calculator
        %hError = comm.ErrorRate;

        %recordar hacer sigma/2 en el teorico tambien
        simbtx = randi(M,nsimbs*length(SNRbdB),1)-1;% simbolos a transmitir (en 1:M);
        
        if modindex~=4
            simbmod = step(modulator,simbtx); % simbolos modulados 
            %simbmod_ruido = step(hAWGN,reshape(simbmod,nsimbs,length(SNRbdB)));%en cada columba una canal distinto
            sigma = sqrt(Eavg./SNR); % desvio de ruido
            ruido = complex(randn(nsimbs,1),randn(nsimbs,1))*sigma;
            simbmod_ruido = simbmod + reshape(ruido,nsimbs*length(SNRbdB),1);

            simbrx= step(demodulator,reshape(simbmod_ruido,nsimbs*length(SNRbdB),1)) ; % simbolos recibidos
            %errores = step(hError, simbtx, simbrx);
        else %FSK -- modulacion ortogonal
            simbmod = zeros(length(simbtx)*M,1); % en grupos de M las coordenadas m-dimensionales
            simbmod([0:nsimbs*length(SNRbdB)-1]'*M+simbtx+1)=1;
            %
            sigma = sqrt(Eavg./SNR); % desvio de ruido
            ruido = complex(randn(nsimbs*M,1),randn(nsimbs*M,1))*sigma;% nM x length(sigma)
            
            simbmod_ruido = reshape(simbmod,nsimbs*M,length(sigma)) + ruido;

            [modulos simbrx]=max(abs(eye(M)*reshape(simbmod_ruido,M,[]))); % simbolos recibidos
            simbrx = simbrx'-1;
        end
        %cada columna es un canal con distinta snr y cada elemento es una señal

        Pes = ones(1,nsimbs)*reshape(simbtx ~= simbrx,nsimbs,length(SNRbdB))/nsimbs;

        % Ahora analizo los errores de bit (el Gray lo hizo el modulador )
        %
        bitstx = de2bi(simbtx,k) ; % convierto los simbolos transmitidos a bits
        bitsrx = de2bi(simbrx,k) ; % convierto los simbolos recibidos a bits
        % comparo los bits y cuento cuantos difieren
        Peb = ones(1,nsimbs)*reshape(xor(bitstx,bitsrx)*ones(k,1),nsimbs,length(SNRbdB)) /k/nsimbs;

        % Graficos
        %
        colours_teo=['--b','--g','--r','--c','--m','--y','--k'];
        colours_montecarlo=['o-b','o-g','o-r','o-c','o-m','o-y','o-k'];
        % Probabidad de error de simbolo
        figure(2*modindex-1)
        semilogy(SNRbdB,Pet(SNR),colours_teo(3*i-2:3*i),'linewidth',1.5);hold on
        semilogy(SNRbdB,Pes,colours_montecarlo(3*i-2:3*i),'linewidth',1.5);
        grid on
        box on
        xlabel('SNR por bit (dB)') ;
        ylabel('Pe')
        legends_pe{2*i-1}=strcat(modulator_name,' Pes teórica');
        legends_pe{2*i}=strcat(modulator_name,' Pes est montecarlo');
        % Probabilidad de error por bit
        figure(2*modindex)
        semilogy(SNRbdB,Pet(SNR)/k,colours_teo(3*i-2:3*i));hold on
        semilogy(SNRbdB,Peb,colours_montecarlo(3*i-2:3*i),'linewidth',1.5);        
        grid on
        box on
        xlabel('SNR por bit (dB)') ;
        ylabel('Pb')
        legends_pb{2*i-1}=strcat(modulator_name,' Peb teórica');
        legends_pb{2*i}=strcat(modulator_name,' Peb est montecarlo');
        %
        %
    end
    %leyendas
    figure(2*modindex-1);hold off;
    legend(legends_pe,'Location','SouthWest');
    %imprimimos
    print_path_pe =sprintf('../Figs/pes_mod_%i.png',modindex);
    print(print_path_pe,'-dpng');
    %leyendas
    figure(2*modindex);hold off;
    legend(legends_pb,'Location','SouthWest');
    %imprimimos
    print_path_pb =sprintf('../Figs/peb_mod_%i.png',modindex);
    print(print_path_pb,'-dpng');
    %
    % limpiamos datos previos
    %clear all;
    if M==16
        figure(cant_modulaciones+cant_constelaciones+1)
        semilogy(SNRbdB,Pet(SNR)/k,colours_teo(3*modindex-2:3*modindex));hold on;
        semilogy(SNRbdB,Peb,colours_montecarlo(3*modindex-2:3*modindex),'linewidth',1.5); 
        grid on
        box on
        xlabel('SNR por bit (dB)') ;
        ylabel('Pb')
        legend_comp_mod{2*modindex-1}=strcat(modulator_name,' Peb teórica');
        legend_comp_mod{2*modindex}=strcat(modulator_name,' Peb est montecarlo');
    end
end
figure(cant_modulaciones+cant_constelaciones+1)
hold off;
legend(legend_comp_mod,'Location','SouthWest');
print('../Figs/peb_comparacion_modulaciones.png','-dpng');
       
